Regresión Logística

Cargamos las librerías que vamos a utilizar

library(tidyverse)
library(tidymodels)
library(modelr)
library(ISLR)
library(GGally)
library(pROC)
library(cowplot)
library(OneR)
library(rlang)
library(caret)
set.seed(1992)

La regresión logística es útil para problemas de predicción de clases. En este caso vamos a tratar de resolver el problema de predecir si una persona va defaultear su deuda de tarjeta de crédito en base a ciertos predictores.

Levantamos el conjunto de datos

Este conjunto de datos proviene de la librería ISLR (Introduction to Statistical Learning Using R) de James, Witten, Hastie y Tibshirani.

# cargamos los datos y observamos su estructura
default <- Default
glimpse(default)

Rows: 10,000
Columns: 4
$ default <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, N...
$ student <fct> No, Yes, No, No, No, Yes, No, Yes, No, No, Yes, Yes, ...
$ balance <dbl> 729.5265, 817.1804, 1073.5492, 529.2506, 785.6559, 91...
$ income  <dbl> 44361.625, 12106.135, 31767.139, 35704.494, 38463.496...

Tiene 4 variables:

• default: La clase que queremos predecir.

• student: Binaria que indica si la persona es estudiante.

• balance: Balance promedio que le queda a la persona luego de sus pagos mensuales.

• income: Ingreso de la persona.

Exploratorias

Analicemos la distribución de la clase.

default %>% 
  group_by(default) %>% 
  summarise(numero_casos=n())

`summarise()` ungrouping output (override with `.groups` argument)

Vemos que estamos trabajando con un problema de clasificación con un claro desbalance de clase.

Realizamos un gráfico exploratorio completo para ver el comportamiento y las relaciones entre las variables. El color rojo designa a quienes no defaultean y el azul a los que sí.

# graficamos con ggpairs coloreando por variable a predecir
g <- ggpairs(default, 1:4,  title = "Correlograma de variables",
        mapping = aes(colour= default)) + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
  theme_bw()
# hacemos un loop para cambiar los colores del gráfico
for(i in 1:g$nrow) {
  for(j in 1:g$ncol){
    g[i,j] <- g[i,j] + 
      scale_fill_brewer(palette="Set1") +  
      scale_color_brewer(palette="Set1")
        }
}
g

¿Qué pueden decir de la relación entre balance y default?

¿Y entre income y balance?

¿Cuáles parecen ser buenas variables para discriminar entre quienes defaultean y quienes no?

Limpieza

Para modelizar va a ser necesario tener la variable default como numérica. Definimos la variable como {0,1} para los valores {“No”,“Yes”}

default <- default %>% 
  mutate(default = case_when(default == "No" ~ 0,
                             default == "Yes" ~ 1))
default

Planteo del Problema

Queremos estimar \(P(Default=Yes|X)=P(X)\) para cada individuo y a partir de ello poder definir un punto de corte para predecir quiénes son los que van a entrar en default.

Regresión lineal

En este caso estamos modelando la probabilidad de la siguiente manera:

\(P(X)= \beta_0 + \sum\limits_{j=1}^p \beta_j X\)

Veamos que tan bueno es el modelo lineal para esto, usando balance como predictor.

test_mco <- default %>% 
              lm(formula = default ~ balance, data = .) 
tdy = test_mco %>% tidy() 
tdy

test_mco %>% glance()

Ambos estimadores son significativos y el test de significatividad global del modelo también es significativo.

Veamos un gráfico de nuestro modelo.

Parece tener bastantes problemas para estimar la probabilidad de default de los individuos. Por ejemplo, vemos que hay varios individuos a los cuales les asigna una probabilidad negativa.

Regresión Logística

Para evitar estos problemas, usamos la función logística.

\(P(Y=1|X)= \frac{e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}{1+e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}\)

El lado derecho se llama expit

Esta función acota el resultado entre 0 y 1, lo cual es mucho más adecuado para modelar una probabilidad.

Luego de hacer algunas operaciones, podemos llegar a la expresión:

\(\log {\frac{P(x)}{1-P(x)}}= \beta_0 + \sum\limits_{j=1}^p \beta_j X\)

El lado izquierdo es el logaritmo de los odds y se llama logit.

Partición del dataset en train y test

Realizamos una partición una partición entre dataset de entrenamiento (70%) y testeo (30%) usando la función initial_split del paquete rsample de tidymodels.

# fijamos semilla
set.seed(44)
# Partición Train y Test, indicando proporción
train_test <- initial_split(default, prop = 0.7)
# armamos dataframe de testeo y entrenamiento
default <- training(train_test)
test <- testing(train_test)
# vemos el contenido
default %>%
  dim_desc() # 7000 filas

[1] "[7,000 x 4]"

test %>%
  dim_desc() # 3000 filas

[1] "[3,000 x 4]"

Creación de fórmulas

Para aplicar la regresión logística primero usamos la función formulas del paquete modelr para crear un objeto que contiene todas las fórmulas que vamos a utilizar. En .response especificamos la variable respuesta de nuestras fórmulas y luego nombramos las fórmulas que queramos armar.

Así, armaremos distintos modelos combinando distintas variables.

logit_formulas <- formulas(.response = ~ default, # único lado derecho de las formulas.
                         bal = ~ balance, 
                         stud = ~ student,  
                         inc = ~ income,  
                         bal_stud = ~ balance + student, 
                         bal_inc = ~ balance + income, 
                         stud_inc = ~ student + income,  
                         full = ~ balance + income + student  
                         )
logit_formulas # observamos el objeto formulas

$bal
default ~ balance

$stud
default ~ student

$inc
default ~ income

$bal_stud
default ~ balance + student

$bal_inc
default ~ balance + income

$stud_inc
default ~ student + income

$full
default ~ balance + income + student

Creación de modelos

Procedemos a crear los modelos a partir de estas fórmulas.

models <- data_frame(logit_formulas) %>% # dataframe a partir del objeto formulas
  mutate(models = names(logit_formulas), # columna con los nombres de las formulas
         expression = paste(logit_formulas), # columna con las expresiones de las formulas
         mod = map(logit_formulas, ~glm(., family = 'binomial', data = default))) # Que estamos haciendo acá? Que vamos a encontrar en la columna?
models

La funcíón glm() nos permite crear un modelo lineal generalizado (Generalized Linear Model). Al igual que la función lm() toma como argumentos una formula y los datos pero también se debe especificar el argumento family: indicamos la distribución del error y la función link que vamos a utilizar en el modelo.

Algunas familias son:

• Binomial: link=logit

• Poisson: link=log

• Gaussiana: link=identidad

Como estamos trabajando con un fenómeno que suponemos tiene una distribución binomial, así lo especificamos en el parámetro family.

Modelos simples

Probamos los primeros tres modelos, aquellos que tienen un único predictor. Usamos la función tidy para obtener los parámetros estimados para estos tres modelos.

models %>% 
  filter(models %in% c('bal','stud','inc')) %>%
  mutate(tidy = map(mod, tidy))  # Qué realizamos en este paso? Que va a tener esta columna?

Para acceder a los elementos de la nueva columna tidy debemos desanidarla.

models %>% 
  filter(models %in% c('bal','stud','inc')) %>%
  mutate(tidy = map(mod, tidy)) %>%  # Qué realizamos en este paso? Que va a tener esta columna?
  unnest(tidy) %>% 
  mutate(estimate=round(estimate,5), # redondeamos valores para facilitar lectura
         p.value=round(p.value,4))

Observamos que todos los modelos tienen coeficientes significativos excepto el coeficiente asociado al Ingreso (income), el cual tiene un p-valor de 0.17.

Interpretación de los coeficientes

Recordando la ecuación para modelar la probabilidad:

\(P(Y=1|X)= \frac{e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}{1+e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}\)

Se observa que ahora las variables ya no tienen una relación lineal con la probabilidad. En este modelo un coeficiente positivo indica que frente a aumentos de dicha variable la probabilidad aumenta, mientras que un coeficiente negativo nos indica lo contrario. Para nuestro caso:

• El coeficiente de 0.0056 de la variable balance indica que cuando el balance aumenta la probabilidad de default aumenta.

• El coeficiente de 0.40886 de la variable studentYes indica que cuando una persona es estudiante su probabilidad de default aumenta respecto a una persona que no lo es.

• El coeficiente de -0.00001 de la variable income indica que cuando el ingreso aumenta la probabilidad de default disminuye.

Recomendamos leer el capítulo de regresión logística de Interpretable Machine Learning: A Guide for Making Black Box Models Explainable de Molnar Cristoph para una discusión más profunda de la interpretación de un modelo de regresión logística.

Modelo completo

Ahora probamos con un modelo que utiliza las tres predictoras.

models %>% 
  filter(models == "full") %>%
  mutate(tidy = map(mod,tidy)) %>%
  unnest(tidy) %>% 
  mutate(estimate=round(estimate,5),
         p.value=round(p.value,4))

En el modelo con las tres variables los coeficientes de studentYes e income han cambiado de signo, por lo cual su efecto sobre la probabilidad de default es inverso al que observamos en los modelos simples. También hay que agregar ahora que el efecto de la variable se observa al estar controlando por las dos covariables restantes.

Asimismo, se observa que el coeficiente de income no resulta significativo tampoco en este modelo (p-valor>0.05).

Evaluación de todos los modelos

Con map() agregamos la función glance para traer información relevante para la evaluación del modelo. Con unnest() accedemos a dicha información. Por último, agregamos una columna con el porcentaje de deviance explicado por cada modelo y ordenamos el dataset según su valor de deviance.

# Calcular las medidas de evaluación para cada modelo
models <- models %>% 
  mutate(glance = map(mod,glance))
# Obtener las medidas de evaluacion de interes
models %>% 
  unnest(glance) %>%
  # Calculamos la deviance explicada
  mutate(perc_explained_dev = 1-deviance/null.deviance) %>% 
  select(-c(models, df.null, AIC, BIC)) %>% 
  arrange(deviance)

El modelo que utiliza las 3 variables es el que minimiza el deviance, aunque el modelo que excluye a la variable income (poco significativa) está muy cerca. Los 3 últimos modelos reducen muy poco el deviance respecto a la deviance nula.

Gráficos de Evaluación

Realizamos los gráficos para el modelo completo y uno de los modelos con mayor deviance (student+income).

Comenzamos agregando las predicciones con augment con el parámetro type="response". La función augment hereda el argumento type.predict de la función predict.

• Si type.predict = 'link' la predicción es en términos de la función link. En nuestro caso son el logaritmo de las odds, es decir, los valores que toma la expresión logit.

• Si type.predict = 'response' la predicción son las probabilidades de que la observación pertenezca a la clase positiva. En nuestro caso, devuelve la probabilidad de la que persona defaultee.

# Añadir las predicciones
models <- models %>% 
  mutate(pred= map(mod,augment, type.predict = "response"))
#Observaciones con probabilidad más baja
models$pred[1]$bal %>% arrange(.fitted) %>% head(10)

#Observaciones con probabilidad más alta
models$pred[1]$bal %>% arrange(desc(.fitted)) %>% head(10)

Guardamos las predicciones para los modelos mencionados.

# Modelo completo
prediction_full <- models %>% 
  filter(models=="full") %>% 
  unnest(pred, .drop=TRUE)
#Modelo malo
prediction_bad <- models %>% 
  filter(models=="stud_inc") %>% 
  unnest(pred, .drop=TRUE)

Violin plots

# graficamos el modelo completo
violin_full = ggplot(prediction_full, aes(x=default, y=.fitted, group=default, fill=factor(default))) + 
  geom_violin() +
  theme_bw() +
  guides(fill=FALSE) +
  labs(title='Violin plot', subtitle='Modelo completo', y='Predicted probability')
# graficamos el modelo malo
violin_bad = ggplot(prediction_bad, aes(x=default, y=.fitted, group=default, fill=factor(default))) + 
  geom_violin() + 
  theme_bw() +
  guides(fill=FALSE) +
  labs(title='Violin plot', subtitle='Modelo malo', y='Predicted probability')
# mostramos ambos
plot_grid(violin_bad, violin_full)

En los gráficos de violin observamos:

• En el eje de abscisas la clase verdadera: Default o No Default (1 o 0).

• En el eje de ordenadas la probabilidad predicha por nuestro modelo.

• El gráfico nos muestra la distribución de la cantidad de observaciones por su clase real y la probabilidad que le asigna nuestro modelo.

¿Cuál parece ser un punto de corte adecuado para cada modelo?

Gráfico de Hosmer-Lemeshow

Se genera una función para realizar un gráfico de Hosmer-Lemeshow para un dataset. Para ello se fijan los siguientes parámetros:

• dataset: conjunto de datos

• predicted_column: columna con la probabilidad predicha

• class_column: columna con la clase a predecir

• possitive_value: valor de la clase a predecir

• bins: cantidad de grupos del gráfico

• color: color de los puntos

• nudge_x: desplazamiento de la etiqueta en el eje x

• nudge_y: desplazamiento de la etiqueta en el eje y

Hosmer_Lemeshow_plot <- function(dataset, predicted_column, class_column, bins, positive_value, color='forestgreen', nudge_x=0, nudge_y=0.05){
  # Asignar los grupos a las observaciones de acuerdo a la probabilidad predicha
  dataset['group'] <- bin(dataset[predicted_column], nbins = bins, method = 'l', labels=c(1:bins))
  # Contar la cantidad de casos positivos por grupo
  positive_class <- dataset %>% filter(!!sym(class_column)==positive_value) %>% group_by(group) %>% count()
  # Obtener la media de las predicciones por grupo
  HL_df <- dataset %>% group_by(group) %>% summarise(pred=mean(!!sym(predicted_column)), count=n()) %>%
            inner_join(.,positive_class) %>%
            mutate(freq=n/count)
  # Gráfico 
  HM_plot <- ggplot(HL_df, aes(x=pred, y=freq)) + 
    geom_point(aes(size=n), color=color) +
    geom_text(aes(label=n),nudge_y = nudge_y)+
    geom_abline(slope = 1, intercept = 0, linetype='dashed') + 
    theme_bw() +
    labs(title='Hosmer-Lemeshow', size='Casos', x="Probabilidad Predicha", y="Frecuencia observada")
  return(HM_plot)
}

Generamos los gráficos pasandole lo parámetros.

# modelo completo
Hosmer_Lemeshow_plot(prediction_full, '.fitted', 'default', 10, 1) +
  labs(subtitle="Modelo completo")

# modelo malo
Hosmer_Lemeshow_plot(prediction_bad, '.fitted', 'default', 10, 1, color = "firebrick", nudge_y = 0.003) + scale_x_continuous(limits = c(0.02,.06)) + scale_y_continuous(limits = c(.02,.06)) + labs(subtitle="Modelo malo")

En los gráficos de Hosmer-Lemeshow observamos:

• En el eje de abscisas la probabilidad predicha de default.

• En el eje de ordenadas la frecuencia de clase, el cociente entre cantidad de individuos que están en default y el total de individuos.

• La línea punteada designa la igualdad entre probabilidad predicha y frecuencia de clase.

• Los círculos, que se construyen de la siguiente manera:

  • Se dividen a las observaciones en bins en base a la probabilidad predicha
  • Se calcula la frecuencia de clase para cada bin
  • En base a estas dos coordenadas se ubica al círculo en el gráfico
  • El número y tamaño indican la cantidad de observaciones en dicho grupo

Aquellos círculos que se ubiquen por encima de la línea punteada indican que el modelo está subestimando la probabilidad para dichos grupos. Mientras que si los círculos se ubican por debajo el modelo está sobreestimando la probabilidad para dichos grupos.

¿Para qué valores parece existir una sobreestimación de la probabilidad? ¿Para cuáles subestimación?

Curvas ROC

# Calculamos curvas ROC
roc_full <- roc(response=prediction_full$default, predictor=prediction_full$.fitted)

Setting levels: control = 0, case = 1
Setting direction: controls < cases

roc_bad <- roc(response=prediction_bad$default, predictor=prediction_bad$.fitted)

Setting levels: control = 0, case = 1
Setting direction: controls < cases

Graficamos ambas en un mismo plot.

ggroc(list(full=roc_full, bad=roc_bad), size=1) + 
  geom_abline(slope = 1, intercept = 1, linetype='dashed') +
  theme_bw() + 
  labs(title='Curvas ROC', color='Modelo')

print(paste('AUC: Modelo completo', round(roc_full$auc,3)))

[1] "AUC: Modelo completo 0.952"

print(paste('AUC: Modelo malo', round(roc_bad$auc,3)))

[1] "AUC: Modelo malo 0.562"

¿Qué significa cada uno de los ejes?

Punto de corte

Hasta ahora hemos evaluado el modelo de manera general, pero el resultado final del modelo debe consistir en asignar a la persona una clase predicha. En nuestro caso debemos establecer un punto de corte según el cual vamos a separar a las personas en quienes defaultean y quienes no.

Probamos varios puntos de corte y graficamos el accuracy, la sensibilidad o recall, la especificidad y la precisión para cada uno de ellos.

Clases predichas / Clases	Positiva	Negativa
Positiva	True Pos	False Pos
Negativa	False Neg	True Neg

Recordemos que:

\(accuracy = \frac{TP+TN}{TP+FP+FN+TN}\)

\(sensitivity = recall = \frac{TP}{TP+FN}\)

\(specificity = \frac{TN}{TN+FP}\)

\(precision = \frac{TP}{TP+FP}\)

prediction_metrics <- function(cutoff, predictions=prediction_full){
  table <- predictions %>% 
    mutate(predicted_class=if_else(.fitted>cutoff, 1, 0) %>% as.factor(),
           default= factor(default))
  
  confusionMatrix(table(table$predicted_class, table$default), positive = "1") %>%
    tidy() %>%
    select(term, estimate) %>%
    filter(term %in% c('accuracy', 'sensitivity', 'specificity', 'precision')) %>%
    mutate(cutoff=cutoff)
  
}

cutoffs = seq(0.01,0.95,0.01)
logit_pred= map_dfr(cutoffs, prediction_metrics)%>% mutate(term=as.factor(term))

ggplot(logit_pred, aes(cutoff,estimate, group=term, color=term)) + geom_line(size=1) +
  theme_bw() +
  labs(title= 'Accuracy, Sensitivity, Specificity y Precision', subtitle= 'Modelo completo', color="")

¿Qué podemos observar en el gráfico?

¿Podemos definir un buen punto de corte? ¿Cuál sería?

¿Por qué la especificidad tiene ese comportamiento?

Dataset de testing

Seleccionamos el modelo completo, ya que es el que maximizaba el porcentaje de deviance explicada y en base a lo que vimos definimos un punto de corte en 0.25 (pueden probar otros)

sel_cutoff = 0.27
# Creamos el modelo
full_model <- glm(logit_formulas$full, family = 'binomial', data = default)
# Agregamos la predicciones al dataset de testeo
table= augment(x=full_model, newdata=test, type.predict='response') 
# Clasificamos utilizamos el punto de corte
table=table %>% 
  mutate(predicted_class=if_else(.fitted>sel_cutoff, 1, 0) %>% as.factor(), 
         default= factor(default))
# Creamos la matriz de confusión
confusionMatrix(table(table$default, table$predicted_class), positive = "1")

Confusion Matrix and Statistics

   
       0    1
  0 2840   55
  1   51   54
                                         
               Accuracy : 0.9647         
                 95% CI : (0.9574, 0.971)
    No Information Rate : 0.9637         
    P-Value [Acc > NIR] : 0.4092         
                                         
                  Kappa : 0.4864         
                                         
 Mcnemar's Test P-Value : 0.7708         
                                         
            Sensitivity : 0.49541        
            Specificity : 0.98236        
         Pos Pred Value : 0.51429        
         Neg Pred Value : 0.98100        
             Prevalence : 0.03633        
         Detection Rate : 0.01800        
   Detection Prevalence : 0.03500        
      Balanced Accuracy : 0.73889        
                                         
       'Positive' Class : 1              
                                         

Segunda parte

Desbalanceo de la clase

Al explorar el dataset vimos que existía un fuerte desbalance de clase. Sólo el 3% de las observaciones pertenecen a personas que defaultearon. Esto puede tener un efecto en las estimaciones del modelo y su clasificación final.

Existen dos maneras sencillas con las cuales podemos trabajar con una clase desbalanceada:

• Sobre-muestreo (oversampling) de la clase minoritaria

• Sub-muestreo (undersampling) de la clase mayoritaria

La función glm puede tomar como argumento una columna (weigths) de ponderadores para poder hacer esto. Podemos asignar pesos mayores a 1 a la clase minoritaria (oversampling) o menores a 1 a la clase mayoritaria (undersampling). En nuestro problema vamos a realizar un sobresampleo de la clase minoritaria.

# Creamos la columna de ponderadores
default <- default %>% mutate(wt= if_else(default==1,20,1))
# Creamos los modelos con la data 'balanceada'
balanced_models <- data_frame(logit_formulas) %>% # dataframe a partir del objeto formulas
  mutate(models = names(logit_formulas), # columna con los nombres de las formulas
         expression = paste(logit_formulas), # columna con las expresiones de las formulas
         mod = map(logit_formulas, ~glm(.,family = 'binomial', data = default, weights = wt))) #Pasamos la columna wt como ponderadores

Vemos las estimaciones de los parámetros para el modelo completo. ¿Existen cambios?

Ahora veamos la evaluación de los modelos ¿Qué pasó con el porcentaje de deviance explicada? ¿Y con la nula?

Setting levels: control = 0, case = 1
Setting direction: controls < cases
Setting levels: control = 0, case = 1
Setting direction: controls < cases

Violin plots, Curvas ROC y AUCs

Realizamos los gráficos de violin, las curvas ROC y calculamos las AUC.

[1] "AUC Modelo completo: 0.952"
[1] "AUC Modelo malo: 0.562"

¿Dónde se ven los cambios más notorios respecto a nuestros modelos anteriores que no tenían en cuenta el desbalance de la clase?

Punto de corte

Volvemos a realizar las pruebas para varios puntos de corte y graficamos el accuracy, la sensibilidad, la especificidad, el recall y la precision para cada uno de ellos.

¿Qué cambios vemos respecto al gráfico anterior?

Dataset de testing

Probamos en el dataset de testing nuestro modelo balanceado. No es necesario que le creemos pesos al dataset de testeo.

Confusion Matrix and Statistics

   
       0    1
  0 2855   40
  1   54   51
                                          
               Accuracy : 0.9687          
                 95% CI : (0.9618, 0.9746)
    No Information Rate : 0.9697          
    P-Value [Acc > NIR] : 0.6506          
                                          
                  Kappa : 0.5043          
                                          
 Mcnemar's Test P-Value : 0.1800          
                                          
            Sensitivity : 0.56044         
            Specificity : 0.98144         
         Pos Pred Value : 0.48571         
         Neg Pred Value : 0.98618         
             Prevalence : 0.03033         
         Detection Rate : 0.01700         
   Detection Prevalence : 0.03500         
      Balanced Accuracy : 0.77094         
                                          
       'Positive' Class : 1